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Abstract 

The bond graph approach to modelling biochemical networks is extended to allow hierarchical 
construction of complex models from simpler components. This is made possible by represent¬ 
ing the simpler components as thermodynamically open systems exchanging mass and energy via 
ports. A key feature of this approach is that the resultant models are robustly thermodynamically 
compliant: the thermodynamic compliance is not dependent on precise numerical values of pa¬ 
rameters. Moreover, the models are reusable due to the well-defined interface provided by the 
energy ports. 

To extract bond graph model parameters from parameters found in the literature, general and 
compact formulae are developed to relate free-energy constants and equilibrium constants. The 
existence and uniqueness of solutions is considered in terms of fundamental properties of stoi¬ 
chiometric matrices. 

The approach is illustrated by building a hierarchical bond graph model of glycogenolysis in 
skeletal muscle. 


‘Corresponding author 


1 



Contents 


1 Introduction 3 

2 Closed and Open systems 5 

2.1 Simple Example: open and closed systems. 7 

2.2 General case. 8 

2.3 Energy flow. 10 

3 Conversion of kinetic data 11 

3.1 Equilibrium Constants and Free-energy Constants. 12 

3.2 Enzyme Catalysed Reactions. 14 

4 Hierarchical modelling 15 

4.1 Modelling issues . 18 

4.2 Stoichiometric Analysis. 21 

4.3 Simulation. 21 

5 Conclusion 23 

A Supplementary Material 29 

A. 1 Reduced-order equations. 29 

A.2 Conversion of kinetic data . 29 

A.3 Mass-action reactions. 31 

A.4 Relation to the Direct Binding Modular Rate Law of Liebermeister et al. [11] .... 32 

A.5 Relation to the Common Modular Rate Law of Liebermeister et al. [11] . 32 

A.6 Relation to the computational model of Lambeth and Kushmerick [23] . 32 

A.7 TPI . 33 

A.8 Hierarchical modelling. 33 

A.9 Simulation. 35 

A. 10 Virtual Reference Environment. 37 


2 



















1 Introduction 


A central aim of Systems Biology is to develop computational models for simulation and analysis of 
complex biological systems. The scale of this challenge is well illustrated by two examples. Firstly, 
the recent whole-cell model of the pathogenic bacterium mycoplasma genitalium [1] which the authors 
claim accounts for the dynamics over the cell cycle of every gene product. Even for this relatively 
simple organism, the model was developed in the form of many submodels describing different cel¬ 
lular processes such as transcription, translation and metabolism; each submodel required a different 
mathematical description and simulation approach, which were coupled together at simulation time. 
Secondly the Virtual Physiological Fluman project [2, 3] “is a concerted effort to explain how each 
and every component in the body, from the scale of molecules up to organ systems and beyond, works 
as part of the integrated whole” 1 . 

In both examples, a key challenge is to devise modelling strategies and frameworks which not only 
handle complexity but also unite different aspects of cell biology, for example: cellular metabolism, 
signalling and electrophysiology. Such systems involve diverse physicochemical domains such as bio¬ 
chemical, chemoelectric and chemomechanical. There is therefore a need for a modelling approach to 
large multi-domain systems. As discussed below, the bond graph modelling formalism has long been 
used to model multi-domain engineering systems; and recently it has been applied to components of 
biochemical systems [4], The purpose of this paper is extend the modelling of biochemical compo¬ 
nents to the modelling of biochemical systems by developing a hierarchical bond graph representation 
of biochemical systems. 

As detailed in Gawthrop and Crampin [4], thermodynamic aspects of chemical reactions have a 
long history in the Physical Chemistry literature. The roles of chemical potential and Gibb’s free 
energy in biochemical system analysis is described in the textbooks of: Hill [5], Beard and Qian [6] 
and Keener and Sneyd [7]. Indeed, if a computational model of a biochemical system is to be securely 
founded in the underlying chemistry, it must fully comply with the fundamental principles of physical 
chemistry as, for example laid out by Atkins and de Paula [8]. 

Thermodynamic compliance for a set of rate equations has been considered in some detail by 
a number of authors [9-11] allowing this theory to be applied over preexisting models and provid¬ 
ing rules to examine new models and ensure that they are constructed to be compliant. However, 
when rules such as the Haldane constraint are used to constrain a set of numerical parameter values, 
numerical errors may still destroy thermodynamical compliance. The bond graph approach to mod¬ 
elling dynamical systems is well established - a comprehensive account is given in the textbooks of 
Gawthrop and Smith [12], Borutzky [13] and Karnopp et al. [14] and a tutorial introduction for control 
engineers is given by Gawthrop and Bevan [15]. It has been applied to chemical reaction networks 
by Oster et al. [16, 17] and to chemical reactions by Cellier [18], Greifeneder and Cellier [19] and 
Gawthrop and Crampin [4]. Thoma and Atlan [20] discuss “osmosis as chemical reaction through 
a membrane”, LeFevre et al. [21] model cardiac muscle using the bond graph approach and Diaz- 
Zuccarini and Pichardo-Almarza [22] discuss the application of bond graphs to biological system in 
general and muscle models in particular. As illustrated in this paper, the bond graph method provides 
an alternative approach to ensuring thermodynamic compliance: a biochemical reaction network built 
out of bond graph components automatically ensures thermodynamic compliance even if numerical 
values are not correct. In this sense, the bond graph approach is robustly thermodynamically com¬ 
pliant (RTC). Conversely, the discipline imposed by bond graph modelling can alert the modeller to 
errors and inconsistencies that may otherwise have escaped attention. This is illustrated in this paper 
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in the context of revising a well-established model of glycogenolysis [23]. 

Thermodynamics distinguishes between: open systems where energy and matter can cross the 
system boundary; closed systems where energy, but not matter, can cross the system boundary; and 
isolated systems where neither energy nor matter can cross the system boundary [8]. Living organisms 
and their subsystems are open as both matter and energy cross the boundaries, thus, mathematical 
models of living organisms and their subsystems should have explicit formulation as an open sys¬ 
tem. In the bond graph approach used for this paper, open subsystem models explicitly contain ports 
through which matter and energy can flow. 

The need for a systems level approach has been termed “understanding the parts in terms of the 
whole” by Cornish-Bowden et al. [24] and the need for mathematical methods to integrate the range of 
cellular processes has been emphasised by Goncalves et al. [25]. Building such system-level models 
from components is simplified if preexisting validated models can be reused and combined in new 
models. A number of standards have been established to distribute re-usable and consistent models 
of biochemical networks [26], including the mark-up languages SBML [27] and CellML [28], and 
reporting guidelines such as the Minimum information requested in the annotation of biochemical 
models (MIRIAM) [29]. A number of models described with SBML and CellML have been pub¬ 
lished, and many of these are stored within the BioModels Database [30] and the Physiome Reposi¬ 
tory [31], More recently developed, SED-ML [32] now allows “simulated experiments” with SBML 
and CellML models to be described. However, despite the increasing availability of models written 
using standardised formats, it is still a difficult and complex task to combine such models within a 
larger system [33, 34]. 

In the context of a thermodynamically compliant reaction network, model reuse requires that in¬ 
terconnecting models with preexisting thermodynamical compliance produces a model where this 
property is inherited. As discussed in this paper, bond graph models for open systems contain energy 
ports which provide the required interconnections thus allowing the models to be combined to pro¬ 
duce a higher-level model which is robustly thermodynamically compliant. This form of hierarchical 
modelling has been used in the bond graph context by Cellier [35] and makes use of bond graph en¬ 
ergy ports [18, 14]. As discussed in more detail by Gawthrop and Bevan [15] and later within this 
paper, the source sensor (SS) bond graph component, introduced by Gawthrop and Smith [36], pro¬ 
vides a convenient notation for such energy ports. Because the system components are connected by 
energy ports, interactions between the components are two-way, thus the distinction between input 
and output is meaningless and a component’s dynamic properties depend on those components which 
are connected to it. Thus the analogy with active electronic components (e.g. diodes and transistors) 
designed to avoid such two-way interactions is not helpful here. Rather, the analogy is with passive 
electronic components (e.g. resistors and capacitors). 

Gawthrop and Crampin [4] also discuss how bond graph models are parameterised such that pa¬ 
rameters relating to thermodynamic quantities are clearly demarcated from reaction kinetics. How¬ 
ever, this alternative parameterisation leads to an important issue: the parameters required by the bond 
graph model are not the same as the parameters normally used to describe a biochemical reaction net¬ 
work. Thus parameters must be converted from one form to another. This conversion is not trivial for 
two reasons: the given parameters must be checked for thermodynamic compliance for a solution to 
exist, and this solution is not necessarily unique. The problem is solved in this paper and illustrated 
using a specific example from the literature. 

Complex systems may be difficult to comprehend theoretically, and the practical derivation of 
properties is often tedious. These issues can be overcome to some extent by providing general repre¬ 
sentations of, and general formulae relating to arbitrarily complex systems. The stoichiometric matrix 
approach [37, 38] is one such methodology; another approach is provided by van der Schaft et al. 
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[39] who investigate the mathematical structure of balanced chemical reaction networks governed by 
mass action kinetics and provide a powerful mathematical notation. This paper provides a general 
representation of an open biochemical reaction network described using bond graphs as well as gen¬ 
eral formulae based on the stoichiometric matrix [37, 38] approach and on the mathematical structure 
and notation of van der Schaft et al. [39]. As discussed by Gawthrop and Crampin [4, § 3], the use 
of the bond graph concept of causality visually and intuitively ties the mathematical concepts of the 
stoichiometric matrix and its subspaces (representing pathways and conserved moieties) [37, 38] to 
the actual chemistry underlying the system as represented by the biochemical system bond graph. 

One purpose of this paper is to introduce Systems Biologists to the bond graph modelling ap¬ 
proach. To this end, a metabolic pathway model from the literature is re-implemented and reexam¬ 
ined from the bond graph point of view. In particular, Lambeth and Kushmerick [23] consider “a 
computational model for glycogenolysis in skeletal muscle” and provide a detailed thermodynami¬ 
cally consistent model. This model has been further embellished by Vinnakota et al. [40] who also 
provide a detailed comparison with experimental data. Beard [41, Chapter 7] uses this model as an 
example for the simulation of complex biochemical cellular systems, and Mosca et al. [42] uses a 
similar model to examine the interaction of the AKT pathway on metabolism. This paper uses the 
Lambeth and Kushmerick [23] model as an exemplar of how the bond graph approach can be applied 
to build a hierarchical model which is RTC. 

The outline of the paper is as follows. § 2 discusses the modelling of open thermodynamical sys¬ 
tems using bond graphs and the bond graph notion of an energy port ; §2.1 introduces the topic using a 
simple example, § 2.2 gives the general case and § 2.3 focuses on energy flows. Conversion of kinetic 
data into the form required for bond graphs is considered in § 3; in particular 3.1 examines the relation 
between equilibrium constants and free-energy constants and gives a general formula. Existence and 
uniqueness issues are discussed and an alternative derivation of the Wegscheider conditions is given 
and § 3.2 discusses the reformulation of enzyme-catalysed reactions. § 4 develops a hierarchical bond 
graph model based on the well-established model of glycogenolysis by [23] to illustrate the concepts 
developed in this paper, and § 5 concludes the paper with suggested directions for further work. 

A virtual reference environment [43] is available for this paper at https : //sourcef orge . 
net/projects/hbgm/. 

Notation Following van der Schaft et al. [39], we use the convenient notation Exp X to denote the 
vector whose ith element is the exponential of the ith element of X and Ln X to denote the vector 
whose ith element is the natural logarithm of the ith element of X. van der Schaft et al. [39] also 
use the element-wise multiplication operator here, for clarity, we use the equivalent Hadamard 
or Schur product with symbol “o”. The Hadamard or Schur operator o [44, Sec. 7.3] corresponds to 
element-wise multiplication of two matrices which have the same dimensions. This is the same as 
the .* notation within matrix manipulation software. Note that K o X = X o K = (diag K)X = 
(diag X)K, where the notation “diag X” is used for the matrix with elements of X along the diagonal. 

2 Closed and Open systems 

This section discusses the modelling of open thermodynamic systems using bond graphs. In partic¬ 
ular, it is shown how an otherwise closed system becomes an open system by the addition of energy 
ports and the standard bond graph components that represent such ports are introduced. A simple 
motivational example is discussed in § 2.1 and the general case is discussed in § 2.2. 
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(b) Example: external flow v e 
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Figure 1: Closed & open systems, (a) & (b) Simple examples corresponding to reaction schemes (1) 
and (2) with imposed concentrations and flows, respectively. The C:A, C:B, and C:C components 
accumulate the species A, B and C\ and the Re:r1, Re:r2 and Re:r3 represent reactions 1, 2 and 3. 
The SS:[F], SS:[R] and SS:[v_e] components are the energy ports converting a closed system to an 
open one. The system inside the dashed box is a closed system; when SS components (representing 
energy ports) are added the system becomes open, (c) General case. The bond symbols —r correspond 
to vectors of bonds; C, IZe, SS , 0 and 1 symbols correspond to arrays of the associated components; 
and TT components represent the intervening junction structure which transmits energy, comprising 
bonds, junctions and TF components. SS : [//.,,] represents a vector of energy ports imposing external 
chemical potentials; SS : [ve\ represents a vector of energy ports imposing external flows. In partic¬ 
ular, with respect to (a) & (b), C corresponds to the species represented by C:A, C:B, and C:C, IZe 
corresponds to the reactions represented by Re:r1, Re:r2 and Re:r3 and TT : iv/ and TT : N[ 
summarise the connections between reactions and species. With respect to (a), SS : \/i e corresponds 
to SS:[F], SS:[R]; with respect to (b), SS : [ng;] corresponds to SS:[v_e]. TT : n[ and TT : AT 
summarise the connections between reactions and ports, and TT : Nj? between species and ports, in 
(a) & (b). 
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2.1 Simple Example: open and closed systems 

As described in Gawthrop and Crampin [4], the reaction scheme: 

A^B B^C C h A (1) 

has the bond graph representation given within the dashed box of Figures 1(a) and 1(b). This is a 
closed system with respect to the three species A, B and C, thus the system equilibrium (where the 
quantities of A, B and C are constant) corresponds to zero reaction flows. 

In contrast, consider the reaction scheme: 

A^B B^C C + F^A + R (2) 

where the reactants F and R correspond to large external pools with effectively fixed concentrations 
Xf and X r . This is represented in Figure 1(a) where two instances of the effort source and flow 
sensor components, SS:[F] and SS:[R], have been appended to the bond graph. Each SS imposes 
a chemical potential (effort in bond graph parlance) upon the system and inherits the flow on the 
corresponding 1 junction. With the sign convention shown, energy flows into the system through 
SS:[F] and out of the system through SS:[R]. 

Although there is a net flow of energy into the system, in this particular case there is no net flow 
of material. The two 1 junctions are connected by an Re component thus the flows are equal. 

Further details are given within Gawthrop and Crampin [4] as to how the mass-action kinetics of 
this closed system can be directly derived from the bond graph: 


Xa = V 3 - Vi 


X b = vi~ V 2 


X C = V 2 - V 3 (3) 


where: x a , x b and x c arc the amounts of species A, B and C in moles; and v\, v 2 and v 3 arc the molar 
flows associated with reactions 1-3, given by: 

V\ = Ki ( K a x a - K b x b ) v 2 = k 2 (K b x b - K c x c ) v 3 = n 3 (K c x c - K a x a ) (4) 


These linear equations can be written in the notation of § 1 using the stoichiometric matrix N: 


(-10 1 \ 

X = NV V = -ko(N t KoX) where N = I 1-10 

\0 1 - 1 / 

with associated mass and flow vectors, and coefficient matrices: 

(x a \ (vA /«A (K a \ 

X = \x b ] V=\v 2 \ «= U 2 K=\K b \ 

\xj \v 3 J \k 3 J \K c J 


(5) 


( 6 ) 


The open system given by the reaction scheme (2) and the bond graph of Figure 1(a) is given by 
Equations (3) but the equation for v 3 of Equations (4) is replaced by: 


V 3 = K] (KfXfK c X c - K r X r K a Xa) 


(7) 


The simple linear expression (5) is also no longer applicable, but as discussed in § 2.2 it can be 
replaced in general terms by: 


Xi = NiV V = no (Exp N fT Ln K o X — Exp AT T Ln K o A") 


where N i = N r i - n{ (8) 
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and where 
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NT 
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n:= *0 0 lj < 10) 


(9) 


A simulation of this simple example is shown in Figure 2 in order to emphasise the main distinctions 
between open and closed systems. 


2.2 General case 

Figure 1(c) is a generalisation of Figures 1(a) and 1(b), giving a schematic representation of the 
structure for an open system. With reference to Figure 1(c), the five TT components represent the 
five flow transformations: 


f f 
= N- v 


< = Nlv 


)l = N[i 


< = Kv 


v r E = N e v e 


( 11 ) 


and thus the five matrices N- , NT, Ni, iVJ & N E define the stoichiometry of the chemical network. 

As emphasised by Gawthrop and Crampin [4] - bonds, TF components and junctions all transport, 
but do not create or destroy, chemical energy. It follows from this that: 


v 


t A{ = vf 


Hi v T A r i = vl 1 Hi 

Using equations (11), equations (12) become: 

fT 


rT. 


[ j t At = 


He 


v T A r e = vl 1 He 


rT. 

e 


( 12 ) 


v T A[ = v t n( Hi yT Ai = v T N£ T Hi v T A{ = v T N[ He v T A r e = v T N£ r He (13) 

As the equations (13) must be true for all v, it follows that the affinities A can be expressed in terms 
of chemical potential p as: 


a{ = n! 


Hi 


= Nf r Hi Al = N£ 


T 

e He 


A r e = N£ He He = N E T Hi 


(14) 


The fact that energy conservation implies Equation (14) is discussed by Cellier [18]. 

With the notation given in § 1 , the Marcelin formula for reaction flows (used by Gawthrop and 
Crampin [4, Eqn. (2.6)] for the scalar case) may be written as: 


Af 

V = K O (Vq + - v 0 ) + N E V E Where V+ = Exp — 

HI 


and V 0 = Exp (15) 


Similarly, the formula for chemical potential h of Gawthrop and Crampin [4, Eqn. (2.3)] in terms of 
the vector K % of free-energy constants is given by: 


where K, = Exp -£= 
HI 


Hi = RT Ln K t o Xi 


(16) 










Figure 2: Simulation of closed & open systems. The parameters arc arbitrarily chosen as: K a = 
Kb = 1, K c = Kf = 2, K r = 0.25, Xf = X r = 1, k\ = 1, K 2 = 2 & K 3 = 3. (a) As this is a closed 
system, the three reaction flows v\ ... v-£ become zero as time increases, (b) As this is an open system, 
the three reaction flows v\... V 3 do not become zero as the system has an external driver. Flowever, 
as the states become constant, the three flows become equal, (c) Although the closed system flows 
become zero the states do not, as the system has a conserved moiety: the three states have a constant 
sum x a + Xb + x c = 4. (d) The three states tend to different values to those in (c); but the moiety is 
still conserved. 
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where p ( - is the vector of standard chemical potentials for the internal species. It is also convenient 
to re-express the external chemical potentials p e as p e = Zi’TLn K e o X e . Combining Equations 
(15)- (16), the reaction flows V are given by: 

V = kl (U+ - V 0 ~) + V E where U+ = Exp iY^Ln K o X U" = Exp jV rT Ln I< o X (17) 

where the composite stoichiometric matrices NT and N r , the composite state X and the vector of 
free-energy constants K are given by: 



In the context of the simple example shown in Figure 1(a) it can be verified that substituting the 
stoichiometric matrices of Equation (9) into Equation (17) leads to Equations (3) and (4). Equation 
(17) expresses the forward and reverse reaction flows U 0 + and Vq~ in terms of species : amounts X 
and free-energy constants K. It is helpful to reexpress Vf and Vq in terms of quantities related to 
reactions. In particular: 


U+ = K f o X f 
Vq = K r o X r 


where K? = Exp N ^ T Ln K and X? = Exp 1 Ln X 


where K r = Exp N rl Ln K 
In the context of the simple example of Figure 1(a) 

A 0 0 0 0 ^ 


and X r = Exp N rl Ln X 


N f = |0 1 0 0 0 

, 00110 ; 


N rT = 


H) 1 0 0 0> 
0 0 10 0 
.1 0 0 0 1 ; 


(19) 

( 20 ) 


( 21 ) 


hence 


I< f = 



K r = 


I<b 

K c 

J<aK r , 


xf = 



x r = 



( 22 ) 


2.3 Energy flow 


With reference to Figure 1(c): the bonds ( 7 ) and transformers (TT) transmit, but do not create 

or destroy, chemical energy; the C components store, but do not create or destroy, chemical energy; 
the reactions (Tie) dissipate, but do not store, chemical energy; and the ports (SS) transmit chemical 
energy in and out of the system. 

The energy dissipated within the Re components represented by Tie is: 


Pr = ^ ~ = V T A where A = A? — A r (23) 

3 = 1 


Using Equations (14), the affinity Ai can be rewritten as: 

Ai = A{ -A r i = ( N f J - jV[ r ) p,i = —Nf fp = -Nf RT In K.X, (24) 
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t (sec) 

(a) Closed system 



Figure 3: Energy flows for the system illustrated in Figure 1(a) and simulated in Figure 2: Pr is the 
external energy flowing into the system. Pc is the energy flowing into the three C components, and 
Pfl is the energy flowing into (and dissipated by) the three Re components, (a) There is no external 
energy source in the closed system so Pr = 0. Conservation of energy gives Pc + Pr = 0: the 
energy flowing out of the C components is dissipated in the Re components, (b) There is an external 
energy source in the open system with q e = /iTLn K e . Initially, energy flows out of the the open 
system but, in the steady state, the energy into the system is balanced by the energy dissipated in the 
Re components and the energy flow into the C components becomes zero. 


Similarly and A e and A can be rewritten as: 

A e = -Njn e = -NjRT In K e X e and A = -iV T /x = -N t RT In KX (25) 

Hence naming the external energy flow into the system as Pr and the energy flow into the 
C components as Pc it follows that: 

P E = RTV T Nj Ln K e o X e (26) 

P c = RTV t nJ Ln Ki o Xi (27) 

and P R = RTV T N T Ln KoX (28) 

These three equations imply energy conservation: 

Pe = Pc + Pr (29) 

Figure 3 is based on data corresponding to the simulations of Figure 2 and illustrates these equations 
for both closed (no external energy source) and open systems. 

3 Conversion of kinetic data 

For the bond graph formulation of this paper to be widely applicable, it is necessary to develop a 
framework for converting kinetic data from the conventional form used with enzyme modelling into 
the parameters required by the bond graph formulation. 

As an example of the issues involved, the reaction A B has a single equilibrium constant 
K eq , whereas the bond graph formulation uses the two thermodynamic constants K a and Kf } . In 
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solution. Similarly, the reaction A ;= /i ==± C ;= A of Equation (1) has three equilibrium constants. 
In this case, the solution for K a , I\t, and K c does not exist unless the product of the three equilibrium 
constants is unity (detailed balance). In § 3.1 a derivation of the formulae that convert equilibrium 
constants to thermodynamic constants for arbitrary networks of reactions is shown. The solution fully 
accounts for the potential of such existence and uniqueness issues. 

Enzyme-catalysed reactions are commonly represented using Michaelis-Menten kinetics [45]. It 
should be noted that when a reaction is assumed to be irreversible, the associated Michaelis-Menten 
kinetics fail to satisfy thermodynamical compliance. For this reason, Michaelis-Menten like kinet¬ 
ics which are thermodynamically compliant have been developed [9-11]. Section 33.2 focuses on 
the enzyme-catalysed reaction formulation of Gawthrop and Crampin [4, §5 (a)] and compares it to 
previously developed methods. 

In the Supplementary Material, this model is shown to be the same as the as the direct binding 
modular rate law of Liebermeister et al. [11] in §A.2A.4, compared with the common modular rate law 
of Liebermeister et al. [11] in §A.2A.5 and compared with the the computational model of Lambeth 
and Kushmerick [23] in §A.2A.6. As well as the equilibrium constant, mass-action reaction kinetics 
also require a rate constant. The corresponding conversion to bond graph form is discussed in § A.2A.3 
of the Supplementary Material. 

3.1 Equilibrium Constants and Free-energy Constants 

With the notation given in § 1, the scalar de Donder formula [46, Equation! 11)] for the ratio Vo °f the 
forward Vq + and backward Vq _ reaction rates of Equation (15) can be written in vector form as: 



A 

or Ln Vq = —— 


(30) 


where A = A* — A r In a similar fashion to (19) & (20) and using Equation (25), Equation (30) can 
be written as: 

V 0 = K v o X v where K v = Exp (-iV T Ln K) and A"” = Exp (-27 T Ln X) (31) 
Substituting Equation (31) into Equation (30) gives the alternative expression for the affinity A: 


A = RT (Ln K v o X v ) = RT (Ln K v + Ln X v ) 


(32) 


In the context of the simple example of Figure 1(a): 



(33) 


As discussed by [7, Eq 1.15], the equilibrium constant of a reaction is given by exp ^4. Com¬ 
bining all reactions, the vector K eq of equilibrium constants is thus: 


K eq = Exp = Exp — 
F RT v RT 


(34) 
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where the subscript 0 indicates the reference state where X is unity. Hence, substituting unit X into 
Equation (31), Equation (34) becomes: 

K eq = K v = Exp (_jv T Ln K) (35) 

Equilibrium corresponds to A = 0 or Vo = l n „ x t hence, using Equation (31), equilibrium also 
corresponds to: 

X v = X eq where K eq X eq = l n „ x l (36) 

Equation (35) gives an explicit expression for the vector K eq , containing ny reaction equilibrium 
constants, in terms of the vector K which contains the free energy constants of nx species. 

However, the transposed stoichiometric matrix N 7 is not normally full rank, and so it is not 
possible to directly use Equation (35) to give I\ in terms of I\ eq . But, as is now demonstrated, the left 
and right null space matrices (as used to detect conserved moieties and flux pathways [37, 38]) lead 
to the solutions (if any) of Equation (35) giving K in terms of K eq . 

The n v x hr right null-space matrix It of N has the property that: 

NR = 0 or R t N t = 0 (37) 

This matrix is used in stoichiometric analysis to analyse metabolic pathways [37, 38]. Here, it is 
reused to examine thermodynamic constraints. Multiplying equation (35) by R T gives: 

R t Ln K eq = 0 (38) 

Equation (38) defines a thermodynamic constraint on the equilibrium constants, and is a form of 
Wegsclieider condition [11, 45]. 

Assuming that the constraint (38) holds, the Moore-Penrose generalised inverse [44, §6.1] of 
— N t can be used to find a solution for K. In particular: 

Ln K 0 = iY f Ln K eq or I < 0 = Exp (iV+Ln K eq ) (39) 

In general, the solution of Equation (39), K = Kq, is not the only value of K satisfying Equation 
(35). The riQ x nx left null-space matrix G (as used to detect conserved moieties [37, 38]) has the 
property that: 

GN = 0 or N t G t = 0 (40) 

Hence a family of solutions of Equation (35) is given by: 

Ln K = Ln Kq + Ln I\\ = Ln KqoK 1 (41) 

where Ln K\ = 6 ,7 Ln k‘\ (42) 

or K\ = Exp (G T Ln k\) (43) 

where k'\ is an arbitrary no x 1 vector. Equation (41) can be rewritten as: 

I< = KqoK\ (44) 

Finally, we show that the affinity A is also unaffected by the choice of k\. From Equations (25) 
and (41)-(44) it follows that: 

A = -N T RTLn (KoX) = RTLn A' 0 o A, oA" = RT(-N t Ln K 0 - N T Ln A, - M T Ln X) 

(45) 

Using (40) & (43), A 7 Ln K\ = 0. Hence A is unaffected by the choice of k\. It follows that the 
energy-based analysis of § 2.3 is also unaffected by the choice of k\. 
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Example The closed system embedded in Figure 1(a) has stoichiometric matrix N, and left and 
right null space matrices G and R given by: 


N = N i = 




G=(l 1 1) 



(46) 


With this value of R, Equation (38) implies that: 


In K{ q + In K eq + In K{ q = 0 or K eq K eq K eq = 1 (47) 

This is the standard “detailed balance” result applied to a three-reaction loop [7, §1.3]. 

Taking the example further, suppose K eq = (1, 0.5, 2) T (which satisfies Equation (47)). Then 
using the GNU Octave implementation pinv to calculate the pseudo inverse for Equation (39) gives 
I\o = (0.79370, 0.79370, 1.58740) T . Noting that nc = 1, k\ is a scalar and may be chosen as 
ki = 1/0.79370 leading to K = (1, 1, 2) T . 

In contrast, the open system of Figure 1(a) has stoichiometric matrix N, and left and right null 
space matrices G and R are given by 


N = 


/-I 
1 
0 
0 
Vo 


0 1 \ 

-1 0 

1 -1 
0 -1 
0 1 / 


A 1 1 0 0\ 
T Vo 0 0 1 l) 


R = 0 (48) 


Because R = 0, the constraint (38) holds for any choice of the three elements of K eq . For example, 
suppose K eq = (1, 0.5, 16) T , then using the GNU Octave implementation pinv for the pseudo 
inverse, Equation (39) gives K 0 = (0.79370, 0.79370, 1.58740, 2.82843, 0.35355) T . Noting that 
no = 2, k\ has two elements and may be chosen, for example, as k\ = (1/0.79370, 2/1.58740 
leading to K = (1, 1, 2, 2, 0.25) T . 


3.2 Enzyme Catalysed Reactions 

The kinetics of enzyme catalysed reactions are usually described by versions of the Michaelis-Menten 
equations; see, for example, Klipp et al. [45, §2.1] or Keener and Sneyd [7, §1.4]. The key issue is 
that the equations should be thermodynamically compliant [11]. 

Gawthrop and Crampin [4, §5(a)] give one possible formulation of the kinetics for enzyme catal¬ 
ysed reactions that is guaranteed to be thermodynamically compliant: 


V = K 


-Ee^O 

l + f*‘ 

rvv 


where a v 


KjU 0 + + k 2 V 0 

Kl + K 2 


(49) 


where V/ + and U 0 - are given by Equation (15), eo is the total enzyme concentration, and K\ and k -2 are 
reaction constants for the reactions relating substrate to complex and complex to product respectively. 
k v is a constant related to the enzyme/complex equilibrium. As is now shown this formulation can be 
rewritten as a reversible Michaelis-Menten equation. Defining: 

Pv = ——— gives a v = (1 - p v )V^ + p v V^ (50) 

+ K ’2 
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With reference to Equation (49), define: 


,. Umax ■ — ,. Vr/iax . _, 

= iim v = -- and v max = urn v = - where v max = Atre c eo 

V+-4-0O 1 — Pv tg - — Pv 


(51) 


then Equation (49) can be rewritten as: 


V? 


- W 


v = 


-71 1 __ _ 'll 

u max iu+ u max u- 

l\, v rvy 

1 + % + ^r 

k v k v 


U 7 + 7 — ^ v 

where kj = -r , k v = — 

(1 pv) Pv 


(52) 


Combining Equations (51) and (52): 


v = 


^o + -^o- 


Vmax (l^" l 7 )) ) 


kv 1 H - 777 ((1 — Pv) + PvVq ) kv + (1 — p v ) + PvVq 


(53) 


Remarks 


1. Equation (53) explicitly shows that this particular form of the reversible Michaelis-Menten equa¬ 
tion has three parameters: v ma x, k v and p v . 

2. When p v = 0, Equation (53) becomes the irreversible Michaelis-Menten equation but with the 
addition of the Vq term to give reversibility. 

3. Using Equation (51), the two parameters v max and p v can be computed from v+ w r and v^ nax as: 


Umax | | _ 

_ v-max _ v max j _ v max v max 

Pv — + — _|_ _ CUIU u max — _|_ _ 

2 _|_ v max Vmax i Vmax Vmax “r Vmax 

Vmax 


(54) 


4. Using Equation (52), the third parameter k v can be computed from: 

kv — (1 Pv) k ^ — p v k v 


(55) 


5. Note that the Haldane equation K eq = Vr ^ ax ^ v j s automatically satisfied by Equation (53). 

Vmaxky 

6. In the limit as k v —> oo, Equation (53) becomes the mass-action flow [4, Equation (2.6) ]: 

v = k (U 0 + — u 0 ") where k = - " Lai 

rZ v 


(56) 


4 Hierarchical modelling 

The bond graph representation for open systems detailed in §2 Figure 1, with the SS ports to connect 
systems, provides a basis to construct hierarchical models of biochemical systems that are robustly 
thermodynamically compliant. This approach is illustrated using a well-established model from the 
literature: “A Computational Model for Glycogenolysis in Skeletal Muscle” presented by Lambeth 
and Kushmerick [23]. Although the model has been further embellished by Vinnakota et al. [40] and 
used as an example in the book of Beard [41], we use information and parameters from the original 
model as a basis for the discussion in this paper. 
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Figure 4: The simplified glycolytic pathway as described by Lambeth and Kushmerick [23, Fig¬ 
ure 1], Reactions/enzymes: GP, glycogen phosphorylase; PGLM, phosphoglucomutase; PGI, phos- 
phoglucose isomerase; PFK, phosphofructo kinase; ALD, aldolase; TPI, triose phosphate isomerase; 
GAPDF1, glyceraldehyde 3-phosphate dehydrogenase; PGK, phosphoglycerate kinase; PGM, phos- 
phoglycerate mutase; ENO, enolase; PK, pyruvate kinase; LDH, lactate dehydrogenase; CK, creatine 
kinase; ADK, adenylate kinase; ATPase, ATPases. Metabolites: GLY, glycogen; P,, inorganic phos¬ 
phate; G1P, glucose-1-phosphate; G6P, glucose-6-phosphate; F6P, fructose-6-phosphate; ATP, adeno¬ 
sine triphosphate; ADP, adenosine diphosphate; FBP, fructose-1,6-biphosphate; DHAP, dihydroxy- 
acetone phosphate; GAP, glyceraldehyde 3-phosphate; NAD, oxidised nicotinamide adenine dinu¬ 
cleotide; NADH, reduced NAD; 13BPG, 1,3-bisphosphoglycerate; 3PG, 3-phosphoglycerate; 2PG, 
2-phosphoglycerate; PEP, phosphoenolpyruvic acid; PYR, pyruvate; LAC, lactate; PCr, phosphocre- 
atine; Cr, creatine; AMP, adenosine monophosphate. Enzyme commission (EC) numbers are shown. 
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Figure 5: Hierarchical Bond Graph Model, (a) As discussed in the text, LamKus02 represents the 
top-level model of the system in Figure 4. (b) GLY2LAC represents one of the three submodels in 
(a). (c)&(d) GLY2FBP&FBP2GAP represent two of the three submodels in (b). (e) The reaction 
GP has two parallel reactions GPa and GPb. 
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Figure 4 shows the simplified glycolysis pathway from Lambeth and Kushmerick [23, Figure 
1] using Systems Biology Graphical Notation (SBGN) [47]. There are many ways to subdivide this 
system to create a hierarchical model. Here, we have chosen to divide the system into three conceptual 
modules: 

1. The primary glycolytic reaction chain leading from glycogen ( GLY ) to lactate {LAC) which 
converts adenosine diphosphate ( ADP ) and inorganic phosphate {P) into adenosine triphosphate 
{ATP ) making use of the energy stored in glycogen. This module is a producer of ATP. 

2. The pah - of reactions catalysed by creatine kinase CK and adenylate kinase ADK involving cre¬ 
atine (Cr), phosphocreatine {PCr) and adenosine monophosphate {AMP) as well as ATP and 
ADP. This module is a buffer of ATP. 

3. The reactions catalysed by numerous ATPases {ATPase\ for a more comprehensive description 
see Lambeth and Kushmerick [23]) which convert ATP into ADP and P, and use the released 
energy to perform work. This module is a consumer of ATP. 

Figure 5 gives a bond graph representation of this top-level decomposition where the three modules are 
represented by the compound bond graph components labelled GLY2LAC, CrAMP and the simple 
reaction bond graph component Re:ATPase. The metabolites ATP , ADP and P flow between these 
three modules as illustrated, forming the overall system model LamKus02 . The module GLY2LAC 
is the most complex of Figure 5, and for this reason it is itself hierarchically decomposed in to three 
further modules (Figure 5(b)), represented by the compound bond graph components: GLY2FBP, 
FBP2GAP and GAP2LAC. Bond graph representations for these reactions and species are given in 
Figure 5. The additional component SS:GLYo is discussed in § 4.1. The module corresponding to 
CrAMP is simpler and is shown using the bond graph representation of reactions and species in the 
Supplementary Material, Figure 9. 

There are two approximations made in our implementation of the model described by Lambeth 
and Kushmerick [23]: 

1. the allosteric modulation of reactions GP and PFK by AMP is ignored. 

2. the reaction kinetics are represented by Equation (53). As discussed in §A.2, the kinetics presented 
by Lambeth and Kushmerick [23] are essentially the common modular rate law of Liebermeister 
et al. [11] and, as mentioned in § 3.2, Equation (53) is essentially the direct binding modular rate 
law of Liebermeister et al. [11], our kinetics differ except for first-order reactions. 

Neither assumption conflicts with our aim of illustrating the creation of robustly thermodynamically 
compliant model. As discussed in § 5, future work will look at the more complicated case. 

4.1 Modelling issues 

The bond graph approach to modelling imposes discipline on the modelling process and thus exposes 
errors and inconsistencies that may otherwise have escaped attention. The process of revising the 
well-established model of Lambeth and Kushmerick [23], and its reuse by other authors illustrate 
this process. In particular, four issues were encountered during the development of these bond graph 
models and deserve special attention: the parallel reactions GPa and GPb forming the overall GP re¬ 
action; the modelling of the conversion of glycogen GLY to glucose-1-P G1P; the reaction catalysed 
by ATPase\ and reaction directionality. 

2 As discussed by Gawthrop and Crampin [4], junctions (such as those appearing in Figure 5(b)) with only two impinging 
bonds could be deleted. They are often left in place for clarity or to make further connections if the model is further refined. 
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The parallel reactions GPa and GPb The parallel reactions GPa and GPb are represented in 
bond graph terms in Figure 5(e). Applying the analysis of §3 3.1 to the 22 x 17 stoichiometric matrix 
N corresponding the entire network gives a 17 x 1 right null-space matrix R with an entry of — 1 
corresponding to GPa and 1 corresponding to GPb. Thus the condition R T Ln K eq = 0 of Equation 
(38) gives: 


ln Kj Pa + ln Kj pb — 0 or Kj Pa — I<J pb (57) 

Equation (57) must be satisfied for thermodynamical compliance, corresponding to the notion that 
these alternative forms of the enzyme are catalysing the same chemical reaction. In fact, Lambeth and 
Kushmerick [23, Table 1] have Kq Pcl = Kq PI) = 0.42 such that condition (57) is satisfied and it is 
therefore possible to convert equilibrium constants of this model to the free-energy constants required 
for bond graph modelling. 

It should be noted, however, that in its original form the model is not robustly thermodynamically 
compliant: if equality (57) was violated due to computing or human error, the compliance would 
not be enforced. In this context, it is interesting to note that such an error exists in the published 
literature, with the paper of Mosca et al. [42] re-using the model of Lambeth and Kushmerick [23]. 
In particular, although on p. 17 Mosca et al. [42] correctly use Kq P(X = 0.42, the value on p. 18 
incorrectly gives K £ Q Pb = 16.62 thus their model is not thermodynamically compliant. A glance at 
[23, Table 1] reveals that Mosca et al. [42] have inadvertently copied the equilibrium constant for 
Phosphoglucomutase instead of that for Glycogen Phosphorylase B. In contrast, a model expressed in 
bond graph form could have incorrect parameters; but it would still be thermodynamically compliant. 
This is the advantage of robust thermodynamical compliance. 

The conversion of glycogen ( GLY ) to glucose-l-P (G1P) Lambeth and Kushmerick [23] discuss 
“uncertainty in the kinetic function and substrate concentration describing the glycogen phosphorylase 
reaction” and embed a simplified model of the reaction within the overall model. Flowever (as they 
explicitly state) the resulting model is stoichiometrically and thermodynamically inconsistent. Ligure 
4, and the underlying Ligure 1 of Lambeth and Kushmerick [23] imply a reaction: 

GLY + P^ G1P (58) 

and this is consistent with their equation (p. 821) for the rate of change of GLY: GLY' = — flux^p. 
Flowever, their equation (p. 822) for (luxc;p = Vqp implies the reaction: 

GLY + P^ GLY + G1P (59) 

as they explain, this arises by equating two versions of the glycogen molecule which differ in length 
by the presence of a single monomer: Glycogen n and Glycogen n _ 1 . Biologically, this reflects the 
large, polymeric nature of glycogen such that a monomer/subunit can be cleaved with little discernible 
effect. It is explicitly stated as an assumption by [23] that Glycogen n /Glycogen n _ 1 is unity. Unlike 
reaction (58), reaction (59) implies a zero rate of change of GLY r : GLY r ' = 0. 

Using the bond graph approach it is not possible to simultaneously implement the rate change of 
GLY corresponding to reaction (58) with the reaction flux implied by (59). In short, this is because 
the stoichiometric matrix associated with reactions (58) and (59) are different and thus the assumption 
of §2 2.2 that bonds and junctions transmit, but do not create or destroy, chemical energy would 
be violated. Conceptually, this is equivalent to mass creation, as glycogen is cleaved by glycogen 
phosphorylase but does not change. In practice the effects would be limited over short simulation 
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times, but with a longer term goal of building reusable modular models, it will ultimately lead to 
models lacking thermodynamic compliance. 

This issue is addressed by introducing the external substance GLYo, representing the difference 
between Glycogen^ and Glycogen n _ 1 , into the GP catalysed reactions (58) and (59) of Lambeth and 
Kushmerick [23] using the bond graph component SS:GLYo (Figure 5(c)): 

GLY + GLYo + P^ GLY + G1P (60) 

Using the notation of Lambeth and Kushmerick [23]; reaction (60) implies that GLYo' = — fluxcrp 
and GLY' = 0, thus combining the two incompatible expressions for GLY' in a stoichiometrically 
and thermodynamically consistent way. This illustrates how the bond graph methodology forces the 
modeller to tackle such issues by not allowing a stoichiometrically and thermodynamically inconsis¬ 
tent model to be constructed. Thus this revision replaces a thermodynamically inconsistent submodel 
by a thermodynamically consistent, and therefore physically plausible, submodel. This physically 
plausible submodel acts as a placeholder until a more complex chemically-correct submodel can be 
devised. This illustrates the role of the bond graph framework to allow incremental refinement of 
individual submodels within a hierarchical system. 

The ATPase catalysed reaction Lambeth and Kushmerick [23] state that “The model design fea¬ 
tures stoichiometric constraints, mass balance, and fully reversible thermodynamics as defined by the 
Haldane relation.” However they do break this feature by choosing the reaction catalysed by ATPase 
to be irreversible. Unlike Lambeth and Kushmerick [23], the reaction catalysed by ATPase (repre¬ 
sented by Re:ATPase in Figure 5) is reversible and corresponds to: 

ATP^ADP + P (61) 

Thus the model in this paper is fully thermodynamically reversible and a closed system. Of course 
this reaction is, for practical purposes, one way. But, to follow the systematic bond graph modelling 
procedure, this fact should be reflected by the choice of reaction parameters rather than by violating 
thermodynamic principles. This point is discussed in detail by Cornish-Bowden and Cardenas [48] 
who state that “Our view is that it is always best to use reversible equations in metabolic simulations 
for all processes apart from exit fluxes ...”. 

However, as will be discussed in § 4.3, this closed system may be converted into an open system 
by injecting external flows of: GLYf, GLYr and LAC in such a way as to make the concentrations 
of these three substances constant. 

The simulation is also arranged such that that flows can be optionally disconnected to examine 
system connectivity. This is used in §44.3 to disconnect the reactions Fout and ATPase as appro¬ 
priate. 

Reaction directionality The TPI reaction within Lambeth and Kushmerick [23] provides an inter¬ 
esting example about definitions of reaction “direction”, which must be specified within an associated 
model description. On page 823, Lambeth and Kushmerick [23] say that “The forward direction of 
TPI is defined as producing dihydroxy acetone phosphate”. This is why the directions of the bonds 
impinging on Re:TPI in Figure 5(d) are in the directions shown. Lambeth and Kushmerick [23] also 
state that the equilibrium constant for TPI is K eq = 0.052. However, this is inconsistent with the 
stated directionality and should be replaced by the reciprocal value (see Supplementary Material, § A.7 
for more details). The discipline imposed by the bond graph model avoids ambiguity with regard to 
reaction direction. 
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4.2 Stoichiometric Analysis 


1 

ADP AMP ATP 

2 

Cr PCr 

3 

NAD NADH 

4 

DPG NAD P2G P3G PEP PYR 

5 

ADP 2ATP P PCr DHAP 2FBP GAP 2DPG P2G P3G PEP F6P G1P G6P 

6 

GLY 


Table 1: Conserved Moieties. The left null matrix of the stoichiometric matrix for this system reveals 
six conserved moieties. The simulation takes account of these using Equation (64) in § A.l. 

As discussed in the textbooks of Palsson [37, 38] and in the bond graph context by Gawthrop 
and Crampin [4] the left null matrix of the system stoichiometric matrix gives in formation about 
conserved moieties. With reference to Table 1 the bond graph model has six conserved moities: 

• CM 1-3 are obvious from the equations and pooled metabolite components are assumed to be 
constant by Lambeth and Kushmerick [23]. 

• CM 4 corresponds to the “total oxidised” part of [23, Eqn. (3)]. 

• CM 5 corresponds to [23, Eqn. (2)] of which they say “Correct conservation of mass within the 
model was proven for both open and closed systems by calculating the total [free] phosphate [note 
the absence of AMP] using the following equation”. 

• CM 6 arises as the net flow into GLY is zero and reflects the assumption by Lambeth and Kushm¬ 
erick [23] that Glycogcn r) /Glycogcn () _ | is unity. 

Using the reduced order equations (Supplementary Material, (64) § A.l), these six CMs are automat¬ 
ically taken into account and even numerical errors cannot cause drift in these CMs. Reduced order 
equations are utilised in the simulations presented in the following §. 

4.3 Simulation 

The bond graph model of Figure 5 with reaction kinetics defined by Equation (53) was compiled 
into ordinary differential equations using the bond graph software MTT (model transformation tools) 
[49]. Free energy constants were obtained using the methods of § 3.1 and the kinetic parameters were 
derived as in § 3.2. The reduced order equations (64) § A.l were simulated using the Is ode solver 
within GNU Octave [50] numerical software with a maximum time step of O.Olmin for two cases: 

Closed system. The flows associated with ATPase and Font were constrained to be zero. Initial 
conditions were defined in Lambeth and Kushmerick [23, Table 3], together with the extra equations 
from page 813 describing the redox potential (~R) and total creatine abundance, respectively: 

* NAD = 1000 X Cr + Xp Cr = 40mM (62) 

Xnadh 

The base units used by Lambeth and Kushmerick [23] are mM for concentration and minutes for 
time. These units have been used in the following figures except that flows are converted to M/s and 
concentrations to M for computing the energy flows in Figure 7. 


21 







t (min) 
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Figure 6: Simulation: equilibria. For each reaction, the ratio Vo = K v o X v (31) of the forward to 
backward reaction flows, is plotted, (a) Closed system: each ratio tends to unity: the steady state is a 
thermodynamic equilibrium, (b) Open system: some ratios tend to a non-unity value: the steady state 
is a not a thermodynamic equilibrium. 
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Figure 7: Simulation: energy flows. Three energy flows are plotted: Pr the external energy flow, 
Pr the energy dissipated in all reactions (including ATPase), PATPase the power consumed by the 
processes represented by the ATPase reaction, (a) There is no external energy flow or ATPase 
reaction flow: Pr tends to zero as the energy in the internal species is used up. (b) There is external 
energy flow and ATPase reaction flow: Pr tends to Pr as the energy in the internal species is used 
up. In this case, at steady-state about 90% of the energy is associated with processes represented by 
the ATPase reaction. 
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Figure 6(a) indicates that this steady-state is an equilibrium, providing a useful “thermodynamic 
validation” of the model as discussed by Lambeth and Kushmerick [23]. Figure 7(a) shows energy 
flows, respectively, which become zero after an initial transient as the systems settles into equilibrium. 

Open system. The flows associated with ATPase and Fout were enabled, and the concentration 
of LACo was fixed at a constant small value by adding an appropriate external flow ve- The AT Pase 
coefficient was chosen to be 0.75 corresponding to the “Moderate exercise” column of Lambeth and 
Kushmerick [23, Table 4]. The final equilibrium state of the closed system simulation was used as the 
initial state. In contrast to Figure 6(a), Figure 6(b) indicates that the steady-state of the open system 
is not an equilibrium as the ratio Vo is n °t unity. Figure 7(b) shows energy flows which settle to non¬ 
zero values. In particular, the dissipated power (including AT Pase) p /,• becomes equal to the external 
energy flow Pe . Energy flows associated with AT Pase are indicated separately showing that about 
90% of the energy is associated with processes represented by the AT Pase reaction. The remainder 
is dissipated as heat in the other reactions. 

Further figures appear in the Supplementary material. 


5 Conclusion 

This paper extends the bond graph approach of Gawthrop and Crampin [4] to allow the hierarchical 
modelling of biochemical systems with reusable subsystems and robust thermodynamic compliance. 
This requires two extensions: the modelling of open thermodynamical systems using energy ports and 
the conversion of standard enzymatic rate parameters to the parameters required by the bond graph. 

The reimplimentation of “A Computational Model for Glycogenolysis in Skeletal Muscle”, orig¬ 
inally presented by Lambeth and Kushmerick [23], in a bond graph formulation verifies the utility 
of the bond graph approach. In particular, the discipline imposed by the bond graph reformulation 
focuses on four potential problems with the original model: lack of robustness caused by the separate 
specification of the identical equilibrium constant for two parallel reactions (GPa & G Pby, stoi¬ 
chiometric inconsistency arising from equating two versions of the glycogen molecule; the use of an 
irreversible reaction (ATPase)', and confusion arising from reaction directionality (TPI). A further 
advantage of bond graph modelling illustrated by this example is the automatic generation of stoi¬ 
chiometric matrices and hence conserved moieties. Currently, these must be identified and imposed 
by model developers and this becomes a laborious process with increasing system sizes. 

The model of enzymatic reactions used in this paper is that previously derived by Gawthrop and 
Crampin [4], The relationship between this particular bond graph formulation and some models of 
enzyme kinetics developed by Liebermeister et al. [11] is explored in the paper. Future work will 
examine the development of bond graph representations for a wider range of enzyme models [51]. 
This would facilitate the inclusion of allosteric modulation, such as that exerted by AMP over the 
GPb reaction. 

Metabolic control analysis (MCA) [52] analyses the feedback control behaviour via sensitivity. 
There is a well-established theory of sensitivity bond graphs [53-55], which we will use to give a 
bond graph interpretation of MCA. A number of authors have discussed the role of control theory in 
systems biology [56-60] and there is a well-established theory of control in the context of bond graphs 
[61-63]. Future work will examine feedback control of biochemical networks from the bond graph 
point of view. Metabolic networks are further controlled over longer time scales by gene expression 
modulation of maximum reaction rate. Thus, future work will also examine the modulation of energy 
flows by gene expression as, for example, found in the Warburg effect [64]. 
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An important feature of bond graphs not utilised in this paper is the ability to interconnect different 
physical domains. Future work will examine chemo-electrical and chenro-mechanical transduction as, 
for example, found in the cardiac myocyte. 
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A Supplementary Material 

A.l Reduced-order equations 

This section includes material from Gawthrop and Crampin [4, §3(c)] about using reduced-order equa¬ 
tions for simulation. 

Given the reaction flows V of Equation (17), the rate of change of the internal states is given by: 

Xi = NiV (63) 

As discussed by number of authors the presence of conserved moieties leads to potential numerical 
difficulties with the solution of Equation (63) [65, 66]. Using the notation of Gawthrop and Crampin 
[4, §3(c)], the reduced-order state x and the internal state X, are given by: 

x = L xX NiV Xi = L Xx x + G x Xi{ 0) (64) 

Equation (64) was used to generate all of the simulation figures in this paper. 

A.2 Conversion of kinetic data 


Reaction 

K eq 

Vmaxf 

Vmaxr 


K? 

ADK 

2.210e+00 

8.800e-01 

- 

8.640e-02 

1.225e-01 

CK 

2.330e+02 

5.000e-01 

- 

1.330e+01 

1.499e-01 

ALD 

9.500e-05 

1.040e-01 

- 

5.000e-02 

2.000e+00 

TPI 

1.923e+01 

1.200e+01 

- 

3.200e-01 

6.100e-01 

ENO 

4.900e-01 

1.920e-01 

- 

1.000e-01 

3.700e-01 

Fout 

1.000e+00 

2.000e+02 

- 

1.000e+06 

1.000e+06 

PGM 

1.800e-01 

1.120e+00 

- 

2.000e-01 

1.400e-02 

GAPDH 

8.900e-02 

1.265e+00 

- 

6.525e-05 

2.640e-06 

LDH 

1.620e+04 

1.920e+00 

- 

6.700e-04 

1.443e+01 

PGK 

5.71 le+04 

1.120e+00 

- 

1.600e-05 

4.200e-01 

PK 

1.030e+04 

1.440e+00 

- 

2.400e-02 

7.966e+00 

GPa 

4.200e-01 

2.000e-02 

- 

4.000e+00 

2.700e+00 

GPb 

4.200e-01 

3.000e-02 

- 

2.000e-01 

1.500e+00 

PGI 

4.500e-01 

- 

8.800e-01 

4.800e-01 

1.190e-01 

PGLM 

1.662e+01 

4.800e-01 

- 

6.300e-02 

3.000e-02 

PFK 

2.420e+02 

5.600e-02 

- 

1.440e-02 

1.085e+01 

ATPase 

2.497e+05 

7.500e+02 

- 

1.000e+06 

1.000e+12 


Table 2: Parameters from Lambeth and Kushmerick [23, Tables 1&2]. 


Using the methods of § 3.1, the equilibrium constants quoted by Lambeth and Kushmerick [23, 
Tables 1&2] (Table 2) were converted into the free-energy constants required by the bond graph for¬ 
mulation and are listed in Table 3. 

Using the methods of § 3.2, the reaction constants quoted by Lambeth and Kushmerick [23] (Table 
2) were converted into the reaction constants required by the bond graph formulation of § 3.2 and are 
listed in Table 4. 

Section A.3 looks at mass-action reactions as used for ATPase, § A.4 looks at the relationship 
of the approach in § 3.2 to the direct-binding modular rate law of Liebermeister et al. [11], § A.5 to 


29 










Species 

K 

Mo 

RT 

ADP 

7.677e-01 

-2.643e-01 

AMP 

2.776e-03 

-5.887e+00 

ATP 

4.692e+02 

6.151e+00 

Cr 

6.174e-01 

-4.822e-01 

P 

2.447e-03 

-6.013e+00 

PCr 

1.620e+00 

4.822e-01 

DHAP 

1.038e+00 

3.718e-02 

FBP 

1.968e-03 

-6.231e+00 

GAP 

1.996e+01 

2.994e+00 

DPG 

1.512e+06 

1.423e+01 

LAC 

1.748e-18 

-4.089e+01 

LACo 

1.748e-18 

-4.089e+01 

NAD 

1.659e+03 

7.414e+00 

NADH 

6.026e-04 

-7.414e+00 

P2G 

2.406e-01 

-1.425e+00 

P3G 

4.330e-02 

-3.140e+00 

PEP 

4.910e-01 

-7.114e-01 

PYR 

7.795e-08 

-1.637e+01 

F6P 

7.792e-04 

-7.157e+00 

G1P 

5.827e-03 

-5.145e+00 

G6P 

3.506e-04 

-7.956e+00 

GLY 

1.000e+00 

0.000e+00 


Table 3: Bond graph species parameters used in simulation (See § 3.1 


Reaction 

Vmax 

ky 

Pv 

ADK 

3.439e+02 

4.398e-02 

6.092e-01 

CK 

2.418e-02 

1.863e-01 

1.000e+00 

ALD 

1.040e+02 

9.840e-05 

2.375e-06 

TPI 

1.082e+03 

5.760e-01 

9.098e-01 

ENO 

1.695e+02 

2.124e-02 

1.169e-01 

Fout 

1.000e+05 

8.738e-13 

5.000e-01 

PGM 

3.136e+02 

2.425e-03 

7.200e-01 

GAPDH 

3.953e+02 

1.653e-03 

6.875e-01 

LDH 

1.096e+03 

1.796e-14 

4.292e-01 

PGK 

3.527e+02 

5.847e+00 

6.851e-01 

PK 

4.494e+01 

2.823e-04 

9.688e-01 

GPa 

1.233e+01 

6.035e-03 

3.836e-01 

GPb 

2.841e+01 

4.635e-04 

5.303e-02 

PGI 

5.674e+02 

5.978e-05 

6.448e-01 

PGLM 

1.337e+01 

1.023e-05 

9.721e-01 

PFK 

4.239e+01 

3.985e-03 

2.430e-01 

ATPase 

6.001e+05 

3.755e+08 

1.998e-01 


Table 4: Bond graph reaction parameters used in simulation (See § 3.2 
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the common modular rate law of Liebermeister et al. [11] and § A.6 to the computational model of 
Lambeth and Kushmerick [23] . 


A.3 Mass-action reactions 

The mass-action formulation of chemical equations reveals key issues encountered in converting ki¬ 
netic data from enzymatic models into the form required by a bond graph model. Enzyme catalysed 
reactions are discussed in § 3.2. 

The mass-action formulation presented by Gawthrop and Crampin [4, Equation 2.6] uses the 
Marcelin formulation rewritten here as: 

v = k (V 0 + - F 0 ~) where L+ = e Af/RT and V 0 ~ = e Ar/RT (65) 

In terms of the reaction A ^ B 


V+ = K aXa 

In terms of the reaction A + B 2C 


Vq — B b x b 


Vq — B a X a l\ b X b 


V 0 ~ = (K c x c ) 2 


One standard way of writing the Mass-action rate of A ^ B 


( 66 ) 


(67) 


v = n eq \ x a - 


1 


K, 


-x b 


eq 


similarly, A + B ^ 2C can be rewritten as 

1 


V = K eq ( X a X b ~ 


K, 


eq 


where K eq = (68) 

A 6 

where I< eq = K ^ b (69) 


Define K? as the constant on the substrate side and K r as the constant on the product side. In the 
case of A ^ B 


K f = K a and K r = K b (70) 

and in the case of A + B ^ 267 

K f = K a K b and I< r = K; (71) 

It follows that: 

f I<f 

K eq = Rf k and K eq = — (72) 

As Kf can be computed from K, which in turn can be deduced as discussed in § 3.1, it follows 
that k can be deduced from n eq . 
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A.4 Relation to the Direct Binding Modular Rate Law of Liebermeister et al. [11] 


It is now shown using an example that Equation (53) is of the same form as the direct binding modular 
rate law of Liebermeister et al. [11]. Consider the enzyme catalysed reaction A + B 2 C In this 
case: 


I/ 0 + = K a K b x a x b , V 0 = K 2 c x\ 


(73) 


Equation (53) is of the form: 


v = 


. “f" KqKjjXqXf, 
J max 


— V, 


- K*x\ 
'max ts— 


1 _l_ q I\bX a Xb i K c X c 

+ K+ + K~ 


(74) 


The direct binding modular rate law is given by Equation (4) of Liebermeister et al. [11] and in the 
notation of this paper is: 


L.H~ X a _ h, — [ X c 

^ TSM rsM rv l 


V = U~ 


Kq 1 KM 


1 _1_ x a X b _l_ / Xc 
1 ^ K™ i<m -r [ K M 


Equations (74) and (75) are identical if we set: 


(75) 


Vmax = uk+ 


Vmax 


k m k m 

K + = 

* K a K b 


K~ = 


K, 


M \ 2 


K r 


(76) 


A.5 Relation to the Common Modular Rate Law of Liebermeister et al. [11] 


However, Equation (53) is not the same as the common modular rate law of Liebermeister et al. [11]. 
In the context of the reaction A + B ^ 2 C, the common modular rate law of Liebermeister et al. [11] 
is of the form: 

/ \ 2 


V = U~ 


Xq Xb _ L— I Xc 

K a * \KM 


(77) 


1 'Xa- I 

1 ' K^ 1 ' KM 

a p 


I Xq Xf) 

^ ’W KM 


,JTc_ I t _Xc_ \ 

‘KM + \Km) 


As discussed by Liebermeister et al. [11], the additional denominator terms imply that Equation (77) 
is not the same as Equation (75) but can be considered an approximation to it. However, in the case 
of reaction A ^ B, the common modular and direct binding modular reaction rates are identical. 


A.6 Relation to the computational model of Lambeth and Kushmerick [23] 

The general enzyme catalysed reaction between two species is given by A B and the corresponding 
rate is written by Lambeth and Kushmerick [23] as: 


v = 


V _ V Xb 

v maxj v maxr tsM 

a v b 

1 _|_ x a I Xb 

1 + I<M + K m 


where V m „.xr = V r 


K m k m 

1X b 1X eq 


max f J\M 


(78) 


In this case VX = K a X a and V 0 = K b X b hence Equation (53) becomes: 


v = 


K a X a - K b X b 




(79) 
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Comparing Equations (78) and (79) and using Equation (54) it follows that they are identical if: 


K. 


M 


1 k v 

1 Pv K a 


and 


1 k v 

Pv K b 


( 80 ) 


or: 

K = (1 - p v )K a K f and k v = p v K b K f (81) 

Having deduced p v from the given data using Equation (54), k v can then be deduced from K% 1 using 
K a , as shown in § 3.1. 

Alternatively comparing Equation (52) with Equation (78) gives: 


kt = 


M 


K 

K a 


and k „ = 


K 

K b 


(82) 


Using the expressions for A::+ and k v (52) gives the same result. 


A.7 TPI 

The equilibrium constant is given as K eq = 0.052. However this gives the wrong value of K eqcombo . 
Because the reaction is specified in the “wrong” direction, it is assumed that K eq should be the recip¬ 
rocal of the given value, ie K eq = 19.23. Beard [41] quotes K eq = 19.87; so this alteration seems to 
be correct. 


A.8 Hierarchical modelling 

The bond graphs of the subsystems GAP2LAC and CrAMP appear in Figures 8 and 9. 


-Tig- 


C:P2G 

\ 


-rls- 


C:PEP 

K 


-°k- 


C:LACo 

K 


Figure 8: Submodel: GAP2LAC 

The ODEs, and corresponding flows, automatically generated from the Bond Graph are given by 
the following equations 
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Figure 9: Submodel: CrAMP 


Xadp 2V adk ^pgk ^pk "F ^pfk H - Vatpase 

X a mp = V a dk 

X-atp = Vadk H - Vpgk H” ^p/c ^pfk ^atpase ^c/c 

X C r ~ ^c/c 

Xp — Vgpa ^ppfr “I - Vatpase ^gapdh 

Xpcr = 

Xdhap = ^a/d H” ^£pz 
X fbp = Vpfk ^Zd 
Xgap — ^aZd ^fpZ Vgapdh 

Xdpg = Vpgk H” Vgapdh 

Xlac Vfout “1“ ^Zd/i 

*Z aco — Vfout 
X n ad — Vgapdh “I - Vldh 
X n .adh Vgapdh ^Idh 

Xp2g — ^eno “I - Vpgm 

Xp3g = Vpgk Vpgm 
Xp e p = -|- Veno 

Xpyr = V^,fc Vici/j 
Xf6p = ^Pffi — ^p/fc 
-^glp = ^ppa 4" ^gpft Vpglm 
XgQp — ^P 9 * “I - Vpglm 

Xgly ~ 0 


(83) 
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Vadk — 

V ck = 
Vald — 
Vtpi = 
^eno — 

Vfout = 

Vpgm = 
Vgapdh 

Vldh = 

Vpgk = 

Vgpa = 

Vgpb = 


Vadk 


_u2 v2 I U 
^adp^adp ^ ^ 


amp &atp Xamp X a tp 


Vpgi 


Vpglm 

Vpfk 




atpase 


kadk H” ^adpXadpP 01 ^ k a mpkatpX a mpXatpPadk ^amp^atp-^amp-^-atp 

{^ck{ kadpkpcrXadpXpcr H“ k a tpk C rX a tpX cr )) 


‘X a tpX c 


( k a dpkp Cr X ac ipXp Cr p c i £ k a tpk cr X a ipX cr p c j^ H- k a tpkc 

(Vald( kdhapkgapXdhapXgap kfbpXfbp)) 
ikald “ 1 “ kdhapkgapXdhapXgapPald ~ kfbpXfbpPald “ 1 “ kfbpXfbp) 
{vtiri. kdhapXdhap "i" kgapXgap )) 

( kdhapXdhapptn kgapXgapptn kg a pXg a p + ktn) 

{Veno (kp2gXp2g kpepXp e p)) 


"F 


{keno kp2gXp2gp i 
iVfout (ki 


eno + kp 2 gX p 2 g 


+ kpepXp e pP e 

Xlac ki aco Xi aco )) 


(kfout ki ac Xi ac pf out + ki ac Xi ac + kiaco-^lacoP f out) 

\Vpgm ( kp2gXp2g + kp3gXp3g)) _ 

( kp2g Xp2(j Ppgm kp^gXp^gPpgm “I - kpSgXpSg + kpg rn j 

(' Vgapdh{ kdpg knadhX-dpgXnadh "F kg a pk n adkpX na dXpXg a p )) 

(kdpg knadhXdpgX na dhPgapdh kg a pk na dkpX na dXpXg a pPg a pdh “I - kg a pk na dkpX na dXpXg a p "F kgapdfi ) 

(vidhi, ki ac k na dXi ac X na( i + k nadh kpyr Xnadh Xpy r ) ) 

( klacknad-^-lac-^-nadPldh + kidh ~ knadhkpyrX na dhXpyrPldh “t - k na dh.kpyrX na dhXpy r ) 

(' v pgk ( ~ kadpkdpgX-adpX-dpg "F k a tpkp3gXp3gX a tp)) 

( k a dp kdpg X a dp Xdpg Ppgk k a dp kdpg X a dp Xdp g k a tp kp^g Xp% g X a fp Pp g ^ kpgk) 

(Vpk( k a dpkp e pX ac ipXp e p + k a tpkpyrXpy r X a tp )) 

{k a dpkpepX a( ipXpepppk k a dp kp e p X a dp Xp e p k a tp kpy r Xpy r X a tp Ppj~ kpfc) 

( kglyXglyVgpa( kglpXg\p ~\~ kpXp )) 

( kglpkglyXglpXglyPgpa kglykpXg[y Xp pgp a + kglykpXg[yXp "F kgp a ) 

( kglyXglyVgpb ( kg\pXg\p + kpXp)) 

( kglpkglyXglpXglypgpl) kglykpXg[yXppgp\) + kg[ykpXglyXp ~\~ kgpb) 

(Vpgi( kfQpXfQp + kgQpXgQp )) 

(kfQpXfQpPpgi kgQpXgQpppgi + kgQpXgQp ~\~ kpgi) 

(Vpglmi kglpXglp + kgQpXgQp )) 


( kglpXglpPpglm — kglpXglp kgQpXgQp ppglm kpglm) 

( Vpfk( k a dpkfbpX ac [pX fbp + k a tpkfQpXfQpX a tp) ) 

(kadpkfbpXadpXjbpPpjh k a tpkfQpXfftpXatpppfk + k a ipkfQpXfQpX a ip + kpjk) 

(‘ Vatpase( k a dpkpX a dpXp + k a tpX a tp )) 

( k a dpkpX a dpXpp a ip ase k a tpX a ipp a ip ase + k a tpX at p + k a tpase) 


(84) 


A.9 Simulation 

Further figures corresponding to Section 4.3 appear in Figures 10, 11 and 12. 

Closed system. Figure 11(a) provides another validity check, showing that the conserved moieties 
are constant; because the reduced order equations (§ A.l) were used, this constraint is automatically 
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Figure 10: Simulated concentrations. Evolution of the concentrations for ATP, ADP, P, GLY and 
LAC corresponds to the situation in Lambeth and Kushmerick [23, Figure 2] except that they use unit 
initial states. Following an initial transient, the species concentrations reach steady-state values for 
both the closed and open systems (cf Figure 2). 
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Figure 11: Conserved Moieties. The sum of each conserved moiety remains constant. 
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Figure 12: Simulated mass flows. Both plots show the input mass flow (though SS:GLYo), the output 
mass flow (though Re:Fout) and the ATP flux (though Re:ATPase) (a) In the closed system, the flows 
arc zero, (b) In the open system, the three flows reach a constant steady state. The final value of ATP 
flow is 5.9mM/min; this is close to to value of 6.1mM/min quoted in the “Moderate exercise” column 
of Lambeth and Kushmerick [23, Table 4] 


enforced and thus numerical drift is avoided. Figure 12(a) shows mass flows within the system. 

Open system. Figure 10(b) shows the evolution of simulated concentrations for ATP, ADP, P, 
GLY and LAC which reach new steady state values due to flows induced by the ATPase reaction. 
As with Figure 11(a), Figure 11(b) shows the conserved moieties are constant for the open system. In 
contrast to Figure 12(a), Figure 12(b) shows mass flow which, in this open-system context settle to 
non-zero values. 

A.10 Virtual Reference Environment 

The software required to generate all of the simulation figures shown in this paper from the bond 
graph representation is packaged in the form of a Virtual Reference Environment [43]. It is available 

at https : //sourceforge .net/pro jects/hbgm/. 
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